Extended Abstract for 32 nd AIAA Applied Aerodynamics Conference, 16-20 June 2014 


Aeroelastic Modeling of a Nozzle Startup Transient 

Ten-See Wang 

NASA Marshall Space Flight Center, Huntsville, Alabama, 35812 
Xiang Zhao 1 , 

Alabama A&M University, Huntsville, Alabama, 35762 
Sijun Zhang* 

ESI CFD, Inc., Huntsville, Alabama, 35806 
and 

Yen-Sen Chen s 

Applied Research Laboratory, Hsinchu 30078, Taiwan 

Lateral nozzle forces are known to cause severe structural damage to any new rocket 
engine in development during test. While three-dimensional, transient, turbulent, chemically 
reacting computational fluid dynamics methodology has been demonstrated to capture 
major side load physics with rigid nozzles, hot-fire tests often show nozzle structure 
deformation during major side load events, leading to structural damages if structural 
strengthening measures were not taken. The modeling picture is incomplete without the 
capability to address the two-way responses between the structure and fluid. The objective 
of this study is to develop a tightly coupled aeroelastic modeling algorithm by implementing 
the necessary structural dynamics component into an anchored computational fluid 
dynamics methodology. The computational fluid dynamics component is based on an 
unstructured-grid, pressure-based computational fluid dynamics formulation, while the 
computational structural dynamics component is developed under the framework of modal 
analysis. Transient aeroelastic nozzle startup analyses at sea level were performed, and the 
computed transient nozzle fluid-structure interaction physics presented. 

Aerospace Engineer, ER42, Fluid Dynamics Branch, Propulsion Structure, Thermal, and Fluids Analysis Division, 
Senior Member AIAA. 

1 Associate Professor, Member AIAA. 

* Technical Fellow, Associate Fellow AIAA. 

5 Senior Research Fellow, Associate Fellow AIAA. 


1 



Nomenclature 


Ci,C 2 ,Cs,C f i= turbulence modeling constants, 1.15, 1.9, 0.25, and 0.09. 

C = damping 

C p = heat capacity 

D = diffusivity 

F yZi Fy : F, = integrated force, and component forces in the lateral direction 
/ = frequency 

FI = total enthalpy 

K = thermal conductivity or stiffness 

k = turbulent kinetic energy 

M = mass 

Q = heat flux 

r = Eigen function 

T = temperature 

t = time, s 

u = mean velocities 

V 2 = S u 

x = Cartesian coordinates or nondimensional distance 

Y = physical displacement 

Z = generalized displacement 

a = species mass fraction 

e = turbulent kinetic energy dissipation rate 

[i = viscosity 

ft, = turbulent eddy viscosity (=pC M k 2 /s) 

^ = damping parameter 

Tl = turbulent kinetic energy production 

p = density 

a = turbulence modeling constants, 0.9, 0.9, 0.89, and 1.15 for Eqs. (2), (4-6). 
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t = shear stress 

O = mode shape matrix 

co = chemical species production rate or natural frequency 

Subscripts 

r = radiation 

g = mesh movement 

t = turbulent flow 


I. Introduction 

Nozzle lateral forces during engine startup and shutdown transients, if not properly managed, are known to cause 
severe structural damage to the engine hardware for almost all liquid rocket engines during their initial development 
[1-4]. Transient nozzle side load is therefore considered a high risk item and a critical design issue. For that reason, 
many research efforts [5-26] have been devoted to understanding the side load physics and their impact on the 
magnitude of side loads. For regeneratively-cooled engines such as the Space Shuttle Main Engine (SSME), the 
peak side load generating physics have been identified as the X shock oscillations across the nozzle lip [7]. For film- 
cooled engines such as the Japanese LE-7A engine and the U.S. J-2X engine under development, the major side load 
generating physics have been associated with the jump of the separation line [3, 8]. Other side load physics such as 
the Free-Shock Separation (FSS)-to-Restricted-Shock Separation (RSS) transition have been mentioned as the 
critical physics for the European Vulcain engine [10]. 

In the aforementioned research efforts, computational fluid dynamics (CFD) have been demonstrated as a 
powerful analysis and design tool in computing and understanding the underlying transient side load physics. And 
most of the CFD efforts, have been focused on computing the side load physics with rigid nozzles and not flexible, 
or aeroelastic nozzles. However, during actual hot-firing of a rocket engine, the nozzle wall or the structure of the 
nozzle, often flexes or deforms in response to the lateral aerodynamic forces. The deformation of the nozzle wall 
simultaneously modifies the aerodynamic flowfield and the lateral forces, which in turn affects the nozzle wall 
deformation. This aeroelastic movement of the nozzle wall, which was not considered in the rigid nozzle modeling, 
is indeed one of the important side load physics and needs to be considered. 
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In early 1990’s, Pekkari [21] proposed a simplified model to study the aeroelastic stability of the first bending 
mode of the Vulcain engine. The model consists of an equation of motion describing the displacement of the nozzle 
wall. For calculation of the aerodynamic load on the wall, a simplified pressure distribution, based on a linearized 
supersonic flow theory, was used before the separation point. After the separation point, ambient pressure was used. 
The separation point was established by assuming ratios of separation pressure to ambient pressure as parametric 
studies. The aerodynamic load calculation of the model was later improved by Ostlund [10], where the pressure 
force applied to the nozzle wall, before the separation point, was extracted from three-dimensional Euler 
simulations, while the separation point was estimated through an empirical criterion. These two pioneering 
researches introduced the concept of aeroelastic modeling into nozzle side load studies. In addition, Ostlund [10] 
showed schematically the two asymmetric modes and six buckling modes, which were first visualized in nozzle tests 
of Tuovila and Land [22], 

To follow up on the aeroelastic analysis of Pekkari and Ostlund, it appears that two improvements are in order. 
That is, to refine the simplified wall pressure distribution or the Euler solution on the fluid side by the more accurate 
full CFD solutions, and to go from addressing the first bending mode on the structure side with more detailed 
computational structural dynamics (CSD) formulations that could cover more major deformation modes. These 
improvements were not realized until 2008, when Zhang, et al [25] embarked on probably the first coupled 
CFD/CSD methodology for aeroelastic nozzle study, by coupling among a full Navier-Stokes solver CFD- 
FASTRAN, a CSD code FEM-STRESS, and an interface code MDICE, on a quasi-steady, two-dimensional nozzle. 
In 2013, the same methodology was further demonstrated to show aeroelastic deformation of nozzle wall by Zhao et 
al. [20], on a quasi-steady, three-dimensional analysis of a J-2S nozzle for 0.1818 s of elapsed time. A short while 
earlier in 2012, Blades, et al. [26] developed a similar methodology to that of Zhao, et al. [20] and Zhang, et al. [25], 
by coupling a full Navier-Stokes solver CHEM, a CSD code Abaquas, and an interface code CSE, to demonstrate 
the aeroelastic deformation of a SSME Block II nozzle transient startup for 0.021 s, started at 0.79 s into the startup 
transient. These developments represent advances over earlier simple aeroelastic nozzle analyses, however, it is 
noted that these analyses ignored chemically reacting flows which plays an important role in all transient nozzle 
physics for engine hot-firing at sea level. In addition, none of these analyses computed at times or chamber pressures 
at or close to the occurrence of any major side load events. For example, transient times in the analysis of Blades, et 
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al. ranged from 0.79 s to 0.811 s, which is far from the occurrence of the first major transient nozzle physics for 


SSME startup transient - the FSS-to-RSS transition. 

Engineers at Marshall Space Flight Center and its contractors have been developing a pressure-based, turbulent, 
chemically reacting flow CFD methodology for computing the rigid-body transient nozzle flows. For example, a 
transient SSME nozzle flow separation computation invoking system modeling for inlet flow history was first 
reported in 1992 [24], The analyses that captured and anchored major side load physics and its peak magnitudes for 
regenerative cooled SSME nozzle [7] and film cooled J-2X nozzle [8] were reported in 2009. That anchored CFD 
methodology has been used to support the design and analysis of the current J-2X engine development and selected 
efforts were published in literatures [8, 9, 19, 23], With the recent advances in aeroelastic modeling of nozzle side 
loads [20, 25, 26], and the recognition of the importance of fluid -structure interaction on transient nozzle physics, 
there is a need to improve our design and analysis methodology by considering the aeroelastic modeling. While 
examining the recent aeroelastic analyses [20, 25, 26], it occurred to us that the CFD and CSD computations were 
carried out in different codes and those codes were only connected through an interface code, thus the analyses can 
only be categorized as loosely coupled aeroelastic solutions [20]. In 2013, Wang, et al. [38] implemented the 
necessary CSD formulations directly into the anchored CFD code, and demonstrated the tightly coupled fluid- 
structure interaction algorithms on the transient startup of the SSME Block I engine at sea level, from 2.8 s to 2.876 
s - inside the region of X shock oscillations across the nozzle lip. In this effort, the computations are carried out from 
the command start to the time of full flow in the nozzle, in order to capture all the major physics during the transient 
nozzle startup. 


II. Computational Methodology 

A. Computational Fluid Dynamics 

The computational fluid dynamics (CFD) methodology is based on a multi-dimensional, finite-volume, viscous, 
chemically reacting, unstructured grid, and pressure-based formulation. Time-varying transport equations of 
continuity, species continuity , momentum, total enthalpy, turbulent kinetic energy, and turbulent kinetic energy 
dissipation were solved using a time-marching sub-iteration scheme and are written as: 
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A predictor and corrector solution algorithm was employed to provide coupling of the governing equations. A 
second-order central-difference scheme was employed to discretize the diffusion fluxes and source terms. For the 
convective terms, a second-order upwind total variation diminishing difference scheme was used. To enhance the 
temporal accuracy, a second-order backward difference scheme was employed to discretize the temporal terms. 
Point-implicit method was used to solve the chemical species source terms. Sub-iterations within a time step were 
used for driving the system of second-order time-accurate equations to convergence. Details of the numerical 
algorithm can be found in Refs [27-30]. 

An extended k-s turbulence model [31] was used to describe the turbulence. A modified wall function approach 
was employed to provide wall boundary layer solutions that are less sensitive to the near-wall grid spacing. 
Consequently, the model has combined the advantages of both the integrated-to-the-wall approach and the 
conventional law-of-the-wall approach by incorporating a complete velocity profile and a universal temperature 
profile [32]. A seven-species, nine -reaction detailed mechanism [32] was used to describe the finite-rate, 
hydrogen/oxygen afterburning combustion kinetics. The seven species are H 2 , 0 2 , H 2 0, O, H, OH, and N 2 . The 
thermodynamic properties of the individual species are functions of temperature. The multiphysics pertinent to this 
study have been anchored in earlier efforts, (e.g., SSME axial force and wall heat transfer [27], SSME startup side 
load and dominant shock breathing frequency [7], J-2X startup and shutdown side loads for a nozzlette configuration 
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[8], nozzle film cooling applications [33], and conjugate heat transfer [34], and separated supersonic flows [7, 8, 19, 
23, 24]). 

The Arbitrary Lagrangian-Eulerian formulation given above in Eqs (1) - (6) is implemented in the flow solver in 
a manner that fully satisfies the Discrete Geometric Conservation Laws [35] for each cell in the grid; namely, 
discrete geometric “closure” of the cell at the start and end of each time step, and exact equality between the sum of 
the discrete volumes swept by all moving faces of the cell and the discrete change in the volume of the cell. 
Satisfaction of the second Discrete Geometric Conservation Law ensures that spurious sources and fluxes are not 
generated by the deformation or motion of cells or cell faces, especially for those adjacent to the moving boundaries. 
This in turn avoids loss of accuracy in the fluid dynamic computations and in the coupling across the interface 
between the fluid dynamic and the structural dynamic portions of the problem [20]. 

B. Computational Structural Dynamics 

The structural dynamics response due to fluid flow actions has been analyzed using direct finite- 
element analysis. The aeroelastic equations of motion of solid bodies are given by 


[m]{t}+[c]{t}+M{t}={t} (7) 

where {F} is the displacement vector, [M] is the mass matrix, [C] is the damping matrix, | K] is the 
stiffness matrix, and {F} is the force vector due to the aerodynamic loads and shear stresses. 

The motion Equation (7) of the structure can be solved using modal approach. On the basis of modal 
decomposition of the structure motion with the eigenvectors of the vibration problem, the displacement, 
velocity and acceleration can be transformed [25, 35] to the generalized displacement, velocity and 
acceleration using a transformation matrix, which can be expressed as the following: 

{r}= [®]{z} ; {r}= [<44 ]r}= [®]{z| w 
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Here [®] is the mode shape matrix containing the eigenvectors, orthonormalized with the mass matrix. 


{z}, jzj, and jzj are the generalized displacement, velocity and acceleration vectors, respectively. The 
eigenvectors are orthogonal to both mass and stiffness matrixes and if Rayleigh damping is assumed, it is 
also orthogonal to the damping matrix. Pre-multiplying Equation (8) by [® ] '. we get 

{z}+ [cpf [c][o]{z}+ [cnfMMjz} = [of {F} (9) 

where 


WM®] = 1 


Equation (9) can be written as n individual equations, one for each mode, as follows 


{z i +2% i a) i z i +a)?z i =r i \ 

V = { f } j 


(10) 


Here a)j is the natural frequency for the ith mode $ is the corresponding damping parameter for that 
mode. The solution to Equation (10) can be obtained for each mode using direct integration algorithm. In 
this effort, since the first several mode shapes are usually important, we chose to compute the first four 
mode shapes: ovalization, bending, triangle, and square, as shown in Fig. 1. In our current formulation, a 
structural dynamics software was employed to extract the Eigen modes and Eigen frequencies. 




Fig. 1. The four mode shapes. 

C. Computational Mesh Dynamics 

The geometric and inertial effects of the motions of deforming structures are fed into the flowfield 
through the varying fluid-dynamic boundary conditions at the surface of structures, which account for the 
locations and the velocities of these surfaces. In addition, in order to maintain boundary and grid 
conformity at the fluid-structure interface, the flowfield grid must be deformed to reflect the motion of the 
fluid-structure interfaces. The new coordinates of the fluid mesh nodes are computed by a computational 
mesh dynamics procedure. In the present work, the flowfield grid is deformed at every fluid-structure data 
exchange, which is carried out once at the end of every global time-step, continuously accommodating the 
deformed shape of the structures in the aeroelastic model. The deformation of the fluid-dynamic grid in 
this work is accomplished by using a spring analog approach, which is applicable for any type of meshes. 

The assumption with the spring analog approach is that the mesh nodes are connected like a network 
of springs. By performing a force balance on each of the “spring elements”, an equilibrium balance is 
sought which provides a smooth mesh. If two elements (nodes) are too close, the spring force will repel 
the nodes away from each other. Since each nodal position depends on its neighboring nodes, the 
boundary deformation effects are felt throughout the mesh domain. Therefore, when the geometries of the 
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boundary nodes are changed, the distributed spring system settles the interior nodes into a new 
equilibrium state which can be modeled by 


= 0 ( 11 ) 

7=1 


for all nodes i in the field. In Eq. (11), Si and 8j are the coordinate displacements of node i and its 
neighboring node j, respectively. /V, is the number of neighboring nodes connected to node i, k,-, is the 
spring stiffness for a given edge between i to j which is taken to be inversely proportionally to the length 
of the edge as 


*« = 


1 . 


- x i ) 2 + (?,■ - yj) 2 + (z-i - zjY 


( 12 ) 


D. Transient Startup Sequences 

Transient system-level simulation is a vital part 
of the computational methodology, because it 
provides the time-histories of the inflow properties 
entering the nozzle. Simply put, the ramp rates, or 
histories, of the inlet pressure, fluid temperature 
and species concentrations play an important role 
in determining the type of flow physics, magnitude 
and duration of the flow physics during the 
transient operations. In other words, the time- 
varying inlet flow properties determine the 
residence times of the flow physics inside the 
nozzle. 

The system-level simulation is based on a 



Fig. 2 Computed thruster chamber inlet properties during 
the start-up transient. 
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lumped, control -volume approach to model the rocket engine as a network of components and sub-components. This 
method of transient system-level modeling has been shown to be effective in simulating the low-frequency, transient 
physics associated with the operation of previous and existing rocket engines [SSME, rocket liquid (RL-10), 
integrated power demonstrator, etc.] and therefore, is an important tool in the design and planning of sequencing the 
transient events of rocket engine operation. Figure 2 shows the major inlet flow properties obtained from the system 
model: the time-varying inlet pressure, temperature, and equivalence ratio profiles. These time-varying inlet 
properties were used at the injector faceplate of the thrust chamber for the CFD computation. Two significant 
pressure rise events can be identified in the inlet pressure history of Fig. 2. The first one occurs at 1.5 s due to 
oxygen prime, while the second one occurs at about 2.4 s, caused by the step opening of the oxygen valves in the 
pre-burners. The inlet temperature history shows a sharp jump at 1.5 s, leveling off after 1.75 s, jumps a little bit again at 
2.4s, and increases linearly until around 3.1 s when it reaches the final temperature. The inlet equivalence ratio history 
shows that the thruster environment is fuel rich throughout the start-up transient, especially in the first 1 .5 s, setting up 
the potential for afterburning. That turns out to be the source of the combustion wave, because the pressure jump at 1.5 s 
increases the reaction rate of afterburning, which leads to the generation of the combustion wave. Afterburning plays an 
important part in the subsequent asymmetric flow physics such as the shock transitions and shock breathings across the 
nozzle lip. As mentioned in the beginning of this section, that the route or history between the starting and end points of 
any of the curves in Fig. 2 influences the flow physics intimately, any simplification on any part of the sequence may run 
the risk of missing or degradation of important flow physics. This SSME start transient process involves thermal-fluid 
physics phenomena and safety-based operating practices that are typical of a conventional liquid hydrogen/liquid 
oxygen rocket engine. 


III. Computational Grid Generation 

This section will be completed when the final manuscript is due. 


IV. Boundary and Inlet Conditions 

Since SSME is a first stage engine, fixed freestream boundary conditions were set corresponding to sea level. No- 
slip condition was specified for the solid walls. Time-varying inlet flow boundary conditions were used at the main 
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combustion chamber (MCC) inlet. These inlet flow properties obtained from the engine system simulation include 
the time varying total pressure, temperature and propellant compositions. The time varying propellant composition 
was preprocessed with the Chemical Equilibrium Calculation 
program [36], assuming the propellants were ignited to reach 
equilibrium composition immediately beyond the injector 
faceplate. The larger than unity equivalence ratio throughout the 



5 s ramp period indicates the SSME is operated at fuel rich 
condition and the inlet composition contains mostly steam and 
excess hydrogen, as shown in Fig. 2. At the start command, or 
time zero, the entire flowfield was initialized with quiescent air. 




The presence of air allows the afterburning with the excess fuel 
which contributes to the side force physics. 


Figure 3 Side views of computed nozzle 
shape and deformation contours. 


V. Results and Discussion 

This section will be completed when the final manuscript is due. Figure 3 shows side views of computed nozzle 
shape and deformation contours at four time slices. 


VI. Conclusion 

This section will be completed when the final manuscript is due. 
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